## Run plotFigures.R before running this one!!!

library(ggplot2)
library(tidyr)
library(dplyr)

inPath <- "/Users/meijian/Work/SIMPLE20181102/Input/"
SIMPLE_grid_data <- read_excel("/Users/meijian/Work/SIMPLE20181102/Input/Map/SIMPLE-All_Africa_Countries_List.xlsx")
species_list <- c("africaneggplant", "bambaragroundnut", "cassava", "cocoyam", "cowpea", "fingermillet", "fonio", 
                  "grasspea", "groundnut", "josephscoat", "lablab", "maize", "mungbean", "okra", "pearlmillet", "pigeonpea", 
                  "plantain", "pumpkin","sesame", "sorghum", "soybean", "sweetpotato", "taro", "tef", "tomato", "yams")
#species_list <- "tomato"
speciesNum <- length(species_list)

filterByCROPGRIDS <- function(inPath, species, country_data) {
  library(ncdf4)
  library(raster)
  
  filtered_data <- data.frame()
  area_threshold <- 5
  mask_raster <- raster(xmn=min(country_data$long)-0.25, xmx=max(country_data$long)+0.25,
                        ymn=min(country_data$lat)-0.25, ymx=max(country_data$lat)+0.25,
                        res= 0.5, crs="+proj=longlat +datum=WGS84")
  
  HarvestAreaFile <- paste0(inPath, "Map/CROPGRIDSv1.05_NC_maps/", "CROPGRIDSv1.05_", species, ".nc")
  
  if (file.exists(HarvestAreaFile)){
    nc_data <- nc_open(HarvestAreaFile)
    lon <- ncvar_get(nc_data, "lon")
    lat <- ncvar_get(nc_data, "lat")
    ha <- ncvar_get(nc_data, "harvarea")
    nc_close(nc_data)
    
    rownames(ha) <- lon
    colnames(ha) <- lat
    
    ha_df <- as.data.frame.table(ha)
    colnames(ha_df) <- c('lon','lat','ha')
    ha_df$lon <- as.numeric(as.character(ha_df$lon))
    ha_df$lat <- as.numeric(as.character(ha_df$lat))
    # convert Ocean = -2, Not assessed = -1 to 0
    ha_df$ha[ha_df$ha < 0] <- 0
    
    mask_ha <- rasterize(ha_df[,c('lon','lat')], mask_raster, ha_df[,'ha'], fun=sum, na.rm=T)
    filterTable_data <- data.frame(cbind(xyFromCell(mask_ha, 1:ncell(mask_ha)), round(values(mask_ha),2)))
    colnames(filterTable_data) <- c('long','lat','Area_sum')
    
    filter_data <- subset(filterTable_data, Area_sum >= area_threshold)
    
    filtered_data <- merge(country_data, filter_data,  by = c("long", "lat"), all.x = FALSE, sort=FALSE)
  }
  else {stop(paste0("AutoGenGridInFile.R: Can't find CROPGRIDS data for ",species,"!"))}
  
  return(filtered_data)
}

harvestArea <- list()
for (species in species_list) {harvestArea <- c(harvestArea, list(filterByCROPGRIDS(inPath, species, SIMPLE_grid_data)))}

## baseline weighted average
merged_list_base <- lapply(1:length(harvestArea), function(i) {merge(harvestArea[[i]], df_all_base[[i]], by = c("long", "lat"))})

## calculate weighted average
weighted_average_base <- lapply(merged_list_base, function(df) {
  colSums(as.data.frame(rep(df[11], times=ncol(df[12:31]))) * df[12:31], na.rm = TRUE)/colSums(as.data.frame(rep(df[11], times=ncol(df[12:31]))) * !is.na(df[12:31]), na.rm = TRUE)
})

year_range <- 2000:2019
yield_base <- as.data.frame(year_range)
for (j in 1:speciesNum) {
  yield_base <- cbind(yield_base, as.data.frame(weighted_average_base[[j]]))
}
if (speciesNum>1){
  yield_base <- rbind(yield_base, c("Mean", colMeans(yield_base[,2:ncol(yield_base)], na.rm = TRUE)))
} else {
  yield_base <- rbind(yield_base, c("Mean", mean(yield_base[,2], na.rm = TRUE)))
}
colnames(yield_base) <- c("Year_base",species_list)
write.csv(yield_base, "./weighted_average_figures/26crops_baseline_weighted_average.csv", row.names=FALSE)

## merge all historical and future list with harvest area
## df_all_both comes from plotFigures.R
merged_list <- lapply(df_all_both, function(df) {lapply(1:length(harvestArea), function(i) {
  if (length(df) == 4*speciesNum){
    start_index <- (i-1)*4 + 1
    end_index <- i*4
    sub_df <- df[start_index:end_index]
    lapply(sub_df, function(df2) merge(harvestArea[[i]], df2, by = c("long", "lat")))
  } else {
    start_index1 <- (i-1)*4 + 1
    end_index1 <- i*4
    start_index2 <- (i-1)*4 + 1 + 4*speciesNum
    end_index2 <- i*4 + 4*speciesNum
    sub_df <- df[c(start_index1:end_index1, start_index2:end_index2)]
    lapply(sub_df, function(df2) merge(harvestArea[[i]], df2, by = c("long", "lat")))
  }
})})

## calculate weighted average
weighted_average <- lapply(merged_list, function(df) {lapply(df, function(df1) {lapply(df1, function(df2) {
  colSums(as.data.frame(rep(df2[11], times=ncol(df2[12:41]))) * df2[12:41], na.rm = TRUE)/colSums(as.data.frame(rep(df2[11], times=ncol(df2[12:41]))) * !is.na(df2[12:41]), na.rm = TRUE)
})})})

## generate output file and make plots
year_hist <- as.data.frame(1990:2019)
year_future <- as.data.frame(2035:2064)
yield_lower <- c(500, 400, 650, 300, 1000, 300, 2500, 300, 5000, 1600)
yield_upper <- c(3500, 1000, 950, 800, 3500, 2000, 6500, 600, 11000, 2800)

file_out <- list()
for (j in 1:speciesNum) {
  yield_hist <- cbind(year_hist, as.data.frame(weighted_average[[1]][[j]]), as.data.frame(rowMeans(as.data.frame(weighted_average[[1]][[j]]), na.rm = TRUE)))
  yield_126 <- cbind(year_future, as.data.frame(weighted_average[[2]][[j]][1:4]), as.data.frame(rowMeans(as.data.frame(weighted_average[[2]][[j]][1:4]), na.rm = TRUE)))
  yield_370 <- cbind(year_future, as.data.frame(weighted_average[[2]][[j]][5:8]), as.data.frame(rowMeans(as.data.frame(weighted_average[[2]][[j]][5:8]), na.rm = TRUE)))
  yield_ac_126 <- rowMeans(as.data.frame(weighted_average[[2]][[j]][1:4])) - mean(colMeans(as.data.frame(weighted_average[[1]][[j]])))
  yield_ac_370 <- rowMeans(as.data.frame(weighted_average[[2]][[j]][5:8])) - mean(colMeans(as.data.frame(weighted_average[[1]][[j]])))
  yield_rc_126 <- (rowMeans(as.data.frame(weighted_average[[2]][[j]][1:4])) - mean(colMeans(as.data.frame(weighted_average[[1]][[j]])))) / mean(colMeans(as.data.frame(weighted_average[[1]][[j]]))) * 100
  yield_rc_370 <- (rowMeans(as.data.frame(weighted_average[[2]][[j]][5:8])) - mean(colMeans(as.data.frame(weighted_average[[1]][[j]])))) / mean(colMeans(as.data.frame(weighted_average[[1]][[j]]))) * 100
  
  # generate error bar for bar plots by GCM
  yield_ac_126_gcm <- as.data.frame(weighted_average[[2]][[j]][1:4]) - colMeans(as.data.frame(weighted_average[[1]][[j]]))
  yield_ac_370_gcm <- as.data.frame(weighted_average[[2]][[j]][5:8]) - colMeans(as.data.frame(weighted_average[[1]][[j]]))
  yield_rc_126_gcm <- (as.data.frame(weighted_average[[2]][[j]][1:4]) - colMeans(as.data.frame(weighted_average[[1]][[j]]))) / colMeans(as.data.frame(weighted_average[[1]][[j]])) * 100
  yield_rc_370_gcm <- (as.data.frame(weighted_average[[2]][[j]][5:8]) - colMeans(as.data.frame(weighted_average[[1]][[j]]))) / colMeans(as.data.frame(weighted_average[[1]][[j]])) * 100
  yield_ac_126_error <- cbind(apply(yield_ac_126_gcm, MARGIN = 1, FUN = min), apply(yield_ac_126_gcm, MARGIN = 1, FUN = max))
  yield_ac_370_error <- cbind(apply(yield_ac_370_gcm, MARGIN = 1, FUN = min), apply(yield_ac_370_gcm, MARGIN = 1, FUN = max))
  yield_rc_126_error <- cbind(apply(yield_rc_126_gcm, MARGIN = 1, FUN = min), apply(yield_rc_126_gcm, MARGIN = 1, FUN = max))
  yield_rc_370_error <- cbind(apply(yield_rc_370_gcm, MARGIN = 1, FUN = min), apply(yield_rc_370_gcm, MARGIN = 1, FUN = max))
  
  colnames(yield_hist) <- c("Year_hist","GFDL","IPSL","MPI","MRI","Mean")
  colnames(yield_126) <- c("Year_126","GFDL","IPSL","MPI","MRI","Mean")
  colnames(yield_370) <- c("Year_370","GFDL","IPSL","MPI","MRI","Mean")
  #colnames(yield_ac_126) <- c("Year_126","GFDL","IPSL","MPI","MRI","Mean")
  #colnames(yield_ac_370) <- c("Year_370","GFDL","IPSL","MPI","MRI","Mean")
  
  # 1 in 10 years yield lower and higher bound
  yield_hist_sorted <- as.data.frame(lapply(yield_hist[2:6], sort))
  yield_hist_bounds <- yield_hist_sorted[c(3,28),]
  logic_126_lower <- as.data.frame(sapply(2:6, function(col_num) {yield_126[, col_num] < yield_hist_bounds[1, col_num - 1]}))
  logic_126_higher <- as.data.frame(sapply(2:6, function(col_num) {yield_126[, col_num] > yield_hist_bounds[2, col_num - 1]}))
  logic_370_lower <- as.data.frame(sapply(2:6, function(col_num) {yield_370[, col_num] < yield_hist_bounds[1, col_num - 1]}))
  logic_370_higher <- as.data.frame(sapply(2:6, function(col_num) {yield_370[, col_num] > yield_hist_bounds[2, col_num - 1]}))
  
  logic_126_lower[logic_126_lower == FALSE] <- NA
  yield_126_lower <- yield_126
  yield_126_lower[, 2:6] <- yield_126_lower[, 2:6] * logic_126_lower
  
  logic_126_higher[logic_126_higher == FALSE] <- NA
  yield_126_higher <- yield_126
  yield_126_higher[, 2:6] <- yield_126_higher[, 2:6] * logic_126_higher
  
  logic_370_lower[logic_370_lower == FALSE] <- NA
  yield_370_lower <- yield_370
  yield_370_lower[, 2:6] <- yield_370_lower[, 2:6] * logic_370_lower
  
  logic_370_higher[logic_370_higher == FALSE] <- NA
  yield_370_higher <- yield_370
  yield_370_higher[, 2:6] <- yield_370_higher[, 2:6] * logic_370_higher
  
  ## absolute change and relative change
  yield_ac <- data.frame(
    category = rep(2035:2064, each=2),
    Value = data.frame(Value = c(rbind(yield_ac_126, yield_ac_370))),
    Scenario = rep(c("SSP126", "SSP370"), times = 30)
  )
  yield_rc <- data.frame(
    category = rep(2035:2064, each=2),
    Value = data.frame(Value = c(rbind(yield_rc_126, yield_rc_370))),
    Scenario = rep(c("SSP126", "SSP370"), times = 30)
  )
  yield_ac_error_lower <- data.frame(
    category = rep(2035:2064, each=2),
    Value = data.frame(Value = c(rbind(yield_ac_126_error[,1], yield_ac_370_error[,1]))),
    Scenario = rep(c("SSP126", "SSP370"), times = 30)
  )
  yield_ac_error_upper <- data.frame(
    category = rep(2035:2064, each=2),
    Value = data.frame(Value = c(rbind(yield_ac_126_error[,2], yield_ac_370_error[,2]))),
    Scenario = rep(c("SSP126", "SSP370"), times = 30)
  )
  yield_rc_error_lower <- data.frame(
    category = rep(2035:2064, each=2),
    Value = data.frame(Value = c(rbind(yield_rc_126_error[,1], yield_rc_370_error[,1]))),
    Scenario = rep(c("SSP126", "SSP370"), times = 30)
  )
  yield_rc_error_upper <- data.frame(
    category = rep(2035:2064, each=2),
    Value = data.frame(Value = c(rbind(yield_rc_126_error[,2], yield_rc_370_error[,2]))),
    Scenario = rep(c("SSP126", "SSP370"), times = 30)
  )
  
  ## GCM comparison
  yield_mean_year <- rbind(colMeans(yield_hist[,2:6]), colMeans(yield_126[,2:6]), colMeans(yield_370[,2:6]))
  yield_min_year <- rbind(apply(yield_hist[,2:6], 2, min, na.rm = TRUE), apply(yield_126[,2:6], 2, min, na.rm = TRUE), apply(yield_370[,2:6], 2, min, na.rm = TRUE))
  yield_max_year <- rbind(apply(yield_hist[,2:6], 2, max, na.rm = TRUE),  apply(yield_126[,2:6], 2, max, na.rm = TRUE), apply(yield_370[,2:6], 2, max, na.rm = TRUE))
  
  yield_mean_year_plt <- data.frame(
    category = rep(c("GFDL","IPSL","MPI","MRI","Mean"), each=3),
    Value = data.frame(Value = c(yield_mean_year)),
    Scenario = rep(c("Historical", "SSP126", "SSP370"), times = 5)
  )
  yield_min_year_plt <- data.frame(
    category = rep(c("GFDL","IPSL","MPI","MRI","Mean"), each=3),
    Value = data.frame(Value = c(yield_min_year)),
    Scenario = rep(c("Historical", "SSP126", "SSP370"), times = 5)
  )
  yield_max_year_plt <- data.frame(
    category = rep(c("GFDL","IPSL","MPI","MRI","Mean"), each=3),
    Value = data.frame(Value = c(yield_max_year)),
    Scenario = rep(c("Historical", "SSP126", "SSP370"), times = 5)
  )
  
  #for box/violin plot
  yield_year_violin <- data.frame(
    category = rep(c("GFDL","IPSL","MPI","MRI","Mean"), each=90),
    Value = data.frame(Value = c(as.matrix(rbind(yield_hist[,2:6],yield_126[,2:6],yield_370[,2:6])))),
    Scenario = rep(rep(c("Historical", "SSP126", "SSP370"), each = 30), times = 5)
  )
  
  crop_out <- cbind(yield_hist, yield_126, yield_370)
  file_out <- c(file_out, list(crop_out))
  write.csv(crop_out, paste0("./weighted_average_figures/", species_list[j], '_weighted_average.csv'), row.names=FALSE)

  ## plots
  custom_colors <- c(Mean = "black", GFDL = "blue", IPSL = "green", MPI = "red", MRI = "orange")
  #custom_colors2 <- c(ensemble_mean_126 = "black", SSP126_GFDL_ESM4 = "blue", SSP126_IPSL_CM6A_LR = "green", SSP126_MPI_ESM1_2_HR = "red", SSP126_MRI_ESM2_0 = "orange")
  #custom_colors3 <- c(ensemble_mean_370 = "black", SSP370_GFDL_ESM4 = "blue", SSP370_IPSL_CM6A_LR = "green", SSP370_MPI_ESM1_2_HR = "red", SSP370_MRI_ESM2_0 = "orange")
  
  # Convert the dataframe to long format
  df_long <- yield_hist %>% pivot_longer(cols = -Year_hist, names_to = "GCMs", values_to = "value")
  p1 <- ggplot(df_long, aes(x = Year_hist, y = value, color = GCMs)) +
    geom_line() +
    scale_color_manual(values = custom_colors) +
    labs(title = paste0(species_list[j],"_historical"), x = "Year", y = "Yield (kg/ha)") +
    #ylim(c(yield_lower[j], yield_upper[j])) +  
    theme_minimal() +
    theme(legend.position = "none")
  ggsave(paste0("./weighted_average_figures/", species_list[j], "_weighted_average_hist.png"), p1, width = 7, height = 4, units = "in", dpi = 300) # single plot
  
  df_long <- yield_126 %>% pivot_longer(cols = -Year_126, names_to = "GCMs", values_to = "value")
  df_long_lower <- yield_126_lower %>% pivot_longer(cols = -Year_126, names_to = "GCMs", values_to = "value")
  df_long_higher <- yield_126_higher %>% pivot_longer(cols = -Year_126, names_to = "GCMs", values_to = "value")
  p2 <- ggplot(df_long, aes(x = Year_126, y = value, color = GCMs)) +
    geom_line() +
    geom_point(data=df_long_lower, aes(x=Year_126, y=value, color = GCMs), shape=6, size=3) + # years lower than 1-in-10 low yield in hist
    geom_point(data=df_long_higher, aes(x=Year_126, y=value, color = GCMs), shape=2, size=3) + # years higher than 1-in-10 high yield in hist    
    scale_color_manual(values = custom_colors) +
    labs(title = paste0(species_list[j],"_SSP126"), x = "Year", y = "Yield (kg/ha)") +
    #ylim(c(yield_lower[j], yield_upper[j])) +  
    theme_minimal() +
    theme(legend.position = "none")
  ggsave(paste0("./weighted_average_figures/", species_list[j], "_weighted_average_126.png"), p2, width = 7, height = 4, units = "in", dpi = 300) # single plot
  
  df_long <- yield_370 %>% pivot_longer(cols = -Year_370, names_to = "GCMs", values_to = "value")
  df_long_lower <- yield_370_lower %>% pivot_longer(cols = -Year_370, names_to = "GCMs", values_to = "value")
  df_long_higher <- yield_370_higher %>% pivot_longer(cols = -Year_370, names_to = "GCMs", values_to = "value")
  p3 <- ggplot(df_long, aes(x = Year_370, y = value, color = GCMs)) +
    geom_line() +
    geom_point(data=df_long_lower, aes(x=Year_370, y=value, color = GCMs), shape=6, size=3) + # years lower than 1-in-10 low yield in hist
    geom_point(data=df_long_higher, aes(x=Year_370, y=value, color = GCMs), shape=2, size=3) + # years higher than 1-in-10 high yield in hist    
    scale_color_manual(values = custom_colors) +
    labs(title = paste0(species_list[j],"_SSP370"), x = "Year", y = "Yield (kg/ha)") +
    #ylim(c(yield_lower[j], yield_upper[j])) +  
    theme_minimal() #+
    #theme(legend.position = "none")
  ggsave(paste0("./weighted_average_figures/", species_list[j], "_weighted_average_370.png"), p3, width = 8, height = 4, units = "in", dpi = 300) # single plot
  
  #absolute change
  p4 <- ggplot(yield_ac, aes(x = factor(category), y = Value, fill = Scenario)) +
    geom_bar(stat = "identity", position = "dodge") +
    geom_errorbar(aes(ymin = yield_ac_error_lower$Value, ymax = yield_ac_error_upper$Value), 
                  width = 0.25, position = position_dodge(1)) +
    labs(x = "Year", y = "Absolute Change (kg/ha)") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) # Set the desired angle here
  ggsave(paste0("./weighted_average_figures/", species_list[j], "_weighted_average_ac.png"), p4, width = 8, height = 4, units = "in", dpi = 300) # single plot

  #relative change
  yield_rc_modified <- yield_rc %>% 
    mutate(Value = case_when(Value > 50  ~ 50, Value < -50 ~ -50, TRUE ~ Value)) # make values beyond bounds as bounds for plotting
  p5 <- ggplot(yield_rc_modified, aes(x = factor(category), y = Value, fill = Scenario)) +
    geom_bar(stat = "identity", position = "dodge") +
    geom_errorbar(aes(ymin = yield_rc_error_lower$Value, ymax = yield_rc_error_upper$Value), 
                  width = 0.25, position = position_dodge(1)) +
    labs(x = "Year", y = "Relative Change (%)") +
    theme_minimal() +
    ylim(c(-50, 50)) + 
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) # Set the desired angle here
  ggsave(paste0("./weighted_average_figures/", species_list[j], "_weighted_average_rc.png"), p5, width = 8, height = 4, units = "in", dpi = 300) # single plot
  
  #GCM comparison
  #barplot
  p6 <- ggplot(yield_mean_year_plt, aes(x = factor(category), y = Value, fill = Scenario)) +
    geom_bar(stat = "identity", position = "dodge") +
    geom_errorbar(aes(ymin = yield_min_year_plt$Value, ymax = yield_max_year_plt$Value), 
                  width = 0.25, position = position_dodge(1)) +
    labs(x = "GCMs", y = "Yield (kg/ha)") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 0, hjust = 1)) # Set the desired angle here
  ggsave(paste0("./weighted_average_figures/", species_list[j], "_weighted_average_GCM_compare.png"), p6, width = 8, height = 4, units = "in", dpi = 300) # single plot
  
  ##box/violin plot
  yield_year_violin$Scenario <- as.factor(yield_year_violin$Scenario)
  theme_set(theme_classic() + theme(legend.position = "top"))
  # Initiate a ggplot
  p7 <- ggplot(yield_year_violin, aes(x = category, y = Value)) + 
  # Add mean points +/- SD
  # Use geom = "pointrange" or geom = "crossbar"
    #geom_violin(trim = FALSE) + 
    #stat_summary(
    #  fun.data = "mean_sdl",  fun.args = list(mult = 1), 
    #  geom = "pointrange", color = "black"
    #) + 
  # Combine with box plot to add median and quartiles
  # Change fill color by groups, remove legend
    #geom_violin(aes(fill = category), trim = FALSE) + 
    #geom_boxplot(width = 0.2)+
    #scale_fill_manual(values = c("#00AFBB", "#E7B800", "#FC4E07", "#A04357", "#C0F357"))+
    #theme(legend.position = "none") + 
    geom_violin(aes(color = Scenario), trim = FALSE, position = position_dodge(0.9) ) +
    geom_boxplot(aes(color = Scenario), width = 0.15, position = position_dodge(0.9)) +
    labs(x = "GCMs", y = "Yield (kg/ha)") +
    scale_color_manual(values = c("#E7B800", "#00AFBB", "#FC4E07"))
  ggsave(paste0("./weighted_average_figures/", species_list[j], "_weighted_average_violin.png"), p7, width = 8, height = 4, units = "in", dpi = 300) # single plot
}

weighted_average_hist <- year_hist
weighted_average_126 <- year_future
weighted_average_370 <- year_future
for (species in species_list) {
  weighted_average_data <- read.table(paste0("./weighted_average_figures/", species, "_weighted_average.csv"), header=TRUE, sep=",", stringsAsFactors=FALSE)
  weighted_average_hist <- cbind(weighted_average_hist, weighted_average_data[,6])
  weighted_average_126 <- cbind(weighted_average_126, weighted_average_data[,12])
  weighted_average_370 <- cbind(weighted_average_370, weighted_average_data[,18])
}
weighted_average_hist <- rbind(weighted_average_hist, c("Mean", colMeans(weighted_average_hist[,2:ncol(weighted_average_hist)], na.rm = TRUE)))
weighted_average_126 <- rbind(weighted_average_126, c("Mean", colMeans(weighted_average_126[,2:ncol(weighted_average_126)], na.rm = TRUE)))
weighted_average_370 <- rbind(weighted_average_370, c("Mean", colMeans(weighted_average_370[,2:ncol(weighted_average_370)], na.rm = TRUE)))
weighted_average_126_rc <- (as.numeric(weighted_average_126[nrow(weighted_average_126),2:ncol(weighted_average_126)]) - as.numeric(weighted_average_hist[nrow(weighted_average_hist),2:ncol(weighted_average_hist)])) / as.numeric(weighted_average_hist[nrow(weighted_average_hist),2:ncol(weighted_average_hist)])
weighted_average_370_rc <- (as.numeric(weighted_average_370[nrow(weighted_average_370),2:ncol(weighted_average_370)]) - as.numeric(weighted_average_hist[nrow(weighted_average_hist),2:ncol(weighted_average_hist)])) / as.numeric(weighted_average_hist[nrow(weighted_average_hist),2:ncol(weighted_average_hist)])
weighted_average_rc <- rbind(c(126, weighted_average_126_rc), c(370, weighted_average_370_rc))

colnames(weighted_average_hist) <- c("Year_hist", species_list)
colnames(weighted_average_126) <- c("Year_126", species_list)
colnames(weighted_average_370) <- c("Year_370", species_list)
colnames(weighted_average_rc) <- c("SSP", species_list)
write.csv(weighted_average_hist, "./weighted_average_figures/26crops_historical_weighted_average.csv", row.names=FALSE)
write.csv(weighted_average_126, "./weighted_average_figures/26crops_future126_weighted_average.csv", row.names=FALSE)
write.csv(weighted_average_370, "./weighted_average_figures/26crops_future370_weighted_average.csv", row.names=FALSE)
write.csv(weighted_average_rc, "./weighted_average_figures/26crops_%change_weighted_average.csv", row.names=FALSE)
